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Abstract 

We develop a simple, fast and predictive model of the hierarchical for- 
mation of galaxies which is in quantitative agreement with observa- 
tions. Comparing simulations with observations we place constraints 
on the density of the universe and on the power spectrum of density 
fluctuations. 
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1 Introduction 

Simulations of galaxy formation usually start with a set of gravitating point 
particles with given initial conditions which are then stepped forward in time 
using huge computer resources. Here we consider a complementary approach 
with a density distribution that is continuous instead of discrete. The model 
is presented in Section 2 and compared with observations in Section 3. Con- 
clusions are collected in Section 4. 



2 Model of galaxy formation 

2.1 The hierarchical formation of galaxies 

Let us study linear departures from a critical universe dominated by non- 
relativistic dark matter. [Q] We consider a cube of side L(t) = a c (t)L and 
expand the density in this cube in a Fourier series. For growing modes with 
5 C < 1 we define the linear approximation of the density of the universe as 
follows: 

Plin {x,t,I) = ——{i + 6 c (£,t)} = ?f{l + 6 c (2,t)} 



PcO 

a 



^ <{ 1 + A(O )a c + J2 \ 6 k\ a c ■ ex P 
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where t is the age of the universe, the subscript denotes present day values, 
p' c0 = uj(£l )p c0 , Pco = 3Hq/(8ttG) is the critical density, and a c (t) = (t/t ) 2 ^ 3 . 
ifj: are random phases. Until Section 3.7 we will take the cosmological con- 
stant A = 0. The functions A(f2 ) and uj(Qo) = [2/(3ifo^o)] 2 are given in 
Tables [TJ and |2|. The sum of the Fourier series is over comoving wavevectors 
that satisfy periodic boundary conditions: 

2-7T 

k — — (to + mj + lk) (2) 

J-'O 

where n,m,l = 0, ±1, ±2, ±3 • • • and 1 < [n 2 + m 2 + l 2 ] 1/2 < I. The density 
in Equation ([[]) is due to Fourier components up to wavevectors of magnitude 
ki = 2ttI/L . Note that the critical density p c = Pco/^l is proportional to 
a~ 3 , the density contrast S c grows in proportion to a c , and the "wavelengths" 
A = 27ra c /\k\ of the Fourier components stretch in proportion to a c (we use 
the word "wavelength" even though (j^) does not describe a wave). Our 
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A(fio) 


0.2 
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Table 1: Function A(fi ) for zero cosmological constant. 
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0.2 
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1.000 
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0.679 


1.2 


1.078 


0.5 


0.783 


1.5 


1.189 


0.7 


0.875 


2.0 


1.364 



Table 2: Function u(Qq) for zero cosmological constant. 



convention for Fourier transforms is: 



h = 77 I $c& *o) exp(-ik ■ x)d 3 x 
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S c (x,t ) = y^%exp(ifc 



I" 



;r 



(2*) 1 



Q — * 

<5gexp(ifc • x)d k 



(4) 



where V = Lq. We take the power spectrum of density fluctuations (as 
defined in (P and extrapolated to today in the linear approximation) to be 
of the form 



P(k) = t^jtK 



Aw T 



(27T) 
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1 + r^iy + w 
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(5) 



with w = k/k eq . Note that P{k) ex fc™ at k <C fc eg and P(/c) oc A;™ -4 at 
fc ^> fc eg - The scale invariant Harrison-Zel'dovich spectrum has n — 1. The 
power spectrum P(fc) is normalized such that 



i | <5 c 2 (f , t )rf 3 x = | P(k)d 3 k = J2 l%l 2 - 



(6) 



Equation (|I|) corresponds to the linear approximation valid when 5 C <1. 
Exact solutions can be obtained at peaks of the density fluctuations (assum- 
ing peaks of spherical symmetry and negligible pressure). When S c = 1.06 in 
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Figure 1: The hierarchical formation of galaxies. Dashed lines: Sum of one, 
two and three Fourier components of the density in the linear approximation. 
When the density in the linear approximation reaches the dotted line a galaxy 
forms. Full lines: Halos of several galaxies with density run p oc r~ 2 . As 
larger wavelengths reach the dotted line, larger galaxies form "swallowing" 
galaxies of previous generations. The peculiar displacements of the galaxies 
have not yet been applied. 



the linear approximation (which has already broken down), the exact solu- 
tion is S = 4.55 and the density fluctuation has reached maximum expansion. 
When S c = 5 m = 1.69 in the linear approximation, the exact solution S di- 
verges, the density fluctuation has collapsed and, in our model, a galaxy has 
formed. 

The generation of galaxies on a computer at a given expansion parameter 
a c proceeds as follows (see Figure p. We begin with the largest possible 
galaxy with 1 = 2. We scan the cube searching for maximums of 5 c (x,t). 
If the maximum at x is not "occupied" by a galaxy (see below), and if 
S c {x,t) > S m = 1.69, we generate a new galaxy at x with the characteristics 
listed below. We increase I — > I + 1 and repeat the calculation to form 
galaxies of a smaller generation, and repeat this process up to I = I max 
corresponding to the smallest galaxies that we wish to generate. 

Note that as time goes on, larger and larger distance scales A/ = 2na c /ki 
become non-linear and new galaxies of larger mass form, "swallowing" up 




Figure 2: Scatter plot of y vs. z. 



galaxies of previous generations. Study Figure [l| again. Galaxy formation 
is therefore an ongoing hierarchical process: today's galaxy clusters are the 
seeds of tomorrow's galaxies. 



2.2 Galaxy characteristics 

We assume that the barionic matter discipates energy and falls to the bot- 
tom of the halo potential wells with density run p oc r™ 2 ||. The galaxies 
generated at "generation" / have the following properties regarded as ap- 
proximate descriptions of complex phenomena. Luminous and total 

(luminous + dark) radius: 
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(see Figure [IJ); luminous and total mass: 



n c , 



-My 



lum 



^Hum 

velocity of circular orbits (if spiral): 



M w -7iR 3 Q(t) Pc (t), 



GM 



or 3-dimensional velocity dispersion (if elliptical) 

Vrms = V( v2 ) = V^O" = J 



3GM 



2R 



(7) 



(9) 



(10) 




Figure 3: Scatter plot of M vs. v c to test the Tully-Ficher relation. 



peculiar velocity QTJ : 
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(ii) 



(the sum is over k up to k = kj_i while the peculiar motion given by (" |13| 
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is less than A/2 = ira c /k); and peculiar acceleration: 
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(12) 



What do we mean by "occupied" ? To generate a galaxy of total radius R 
at position a; we require that the distance from x to already generated galaxies 
of radius R' be greater than R' + 0.7 R. The factor ^0.7 was included to 
approximately fill space. 

Finally, after generating all galaxies, we correct their positions x by their 
peculiar motions 
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This step is necessary in order to obtain the galaxy-galaxy correlation. This 
completes the presentation of the model. 

2.3 The simulations 

Our simulations are limited by computer resources (one 500 MHz, 32 bit 
processor). Therefore, we fix Lq = 92Mpc. Even then the contributions to 
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Figure 4: Scatter plot of R vs. v c to test the Samurai relation. 



Pu n from Fourier components with 1=1 and edge effects are non-negligible. 
More realistic simulations require a larger L . We fix /i — 0.6, rj = 2.04 
and k eq = 0.155Mpc _1 ■ VL hl ■ {3.36/N^f} 1 / 2 . This value of k eq corresponds 
to a distance to the horizon equal to ~ 0.5A eg at the time t eq when the 
densities of radiation and matter were equal. N e ff is the effective value 
of Nf, + iNf, where iV 6 (Nf) is the number of boson (fermion) degrees of 
freedom. N e ff = 3.36 for three light neutrino species. The values of k eq and 
rj were obtained from a fit to the "CHDM" model of reference which is in 
agreement with observations. 



3 Comparison with observations 

The results of the model are presented and compared with observations in 
Figures |2| to |^ and Table [|. These figures correspond to the simulation with 
tt = 1.0, n = 1.00 and A = 1706Mpc 3 . Explanations follow. 



3.1 Tully-Fisher, Faber-Jackson and Samurai relations 

The model satisfies M\ um ex v\ oc vf ms , which, if light traces mass, is in 
reasonable agreement with the Tully-Fisher relation for spiral galaxies (L oc 
v% with jjl ~ 4[[| for the infrared luminosity and fi w 2.4 to 2.8 for the blue 
luminosity), or the Faber-Jackson relation for elliptical galaxies (L oc v^ ms 
with \x = 4 ± 10). -^is the galaxy absolute luminosity. The model satisfies 



OC V r 



which is in reasonable agreement with the Samurai relation 
for elliptical galaxies (Ri U m oc vf ms with £ w 1.2 to 1.3). See Figures ^ and[|. 
Elliptical galaxies satisfy the following relation with "remarkably small" H] 
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Figure 5: Histogram of log 10 (M) to test the Schechter distribution. 



scatter: v r . 



oc L 062 R 



-0.50 

lum 



„ A/fO.50 d-0.50 

^ lvI lum ri lum ■ 



This relation is in excellent 



agreement with our model. 

The model satisfies the following useful relations at present: 
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(14) 



(15) 



Both equations are in reasonable agreement with observations if the density 
of matter in optically bright baryons is fli um ~ 0.005. 

3.2 The Schechter distribution 

The observed galaxy luminosity distribution is given by the Schechter rela- 
tion. Q] Assuming that light traces mass, we obtain the following total (lu- 
minous + dark) mass distribution of galaxies: 



dn 



dHM)*** eXP {-M* 



M 



(16) 



where dn is the number of galaxies per unit volume with total mass between 
M and M + dM, 0* = 0.010 exp(±0.4)^Mpc" 3 and M* = Mj um *O O^ = 
2.8-10 13 exp(±0.4)r2o^o~ 1 ^««-[i] The corresponding velocity of circular orbits 
is v c * ■ See Figures [| and ^ 

We choose galaxies of total mass M* to be of "generation" I m 10. 
Then our simulations require L ps I{QM„/t:VLqp c q) 1 ^. We therefore set 
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Figure 6: Histogram of the velocity of circular orbits v c 



L = 92Mpc as indicated above. We also set I max = 30 so that galaxies 
are generated down to a mass ~ M*/27. Galaxies with less mass contribute 
negligibly to the density of the universe and the computer resources needed 
to generate them become prohibitive. The simulations then probe the power 
spectrum in the wavevector range O.lMpc -1 < k < 2.0Mpc _1 . 

3.3 Galaxy-galaxy correlation 

The observed joint probability of finding two galaxies in two volume elements 
separated by r is 



dP oc 



r 



dV 1 dV 2 



(17) 



with 7 = 1.77 ± 0.04 and r = (5.4 ± 1.0)/i " 1 Mpc.[|l Due to the linear 
approximation of x pec we expect agreement with observations to break down 
at small r. See Figure 0. 



3.4 Fluctuation in galaxy counts 

The observed fluctuation in galaxy counts in randomly placed spheres of 
radius r is[|] 
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for r = 9.2h 1 Mpc = 92Mpc/6 which is a convenient radius in view of our 
simulation volume. 
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Figure 7: Histogram of distances between galaxies to test the galaxy-galaxy 
correlation. 



3.5 Fluctuations of the Cosmic Microwave Background 

The Cosmic Microwave Background (CMB) radiation propagates freely since 
matter and radiation decoupled at Td ec ~ 3000K, or at ad ec = (1 + -^ec) -1 ~ 
1100" 1 . The relation between the comoving length at decoupling (or at any 
time with a<l) with the corresponding angle today is[|J 
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Fluctuations on scales 9 < 6' are erased due to the thickness of the last 
scattering surface. The size of the horizon at decoupling, 3ctd ec , corresponds 



1.7°. The size of the horizon at t eg , ~ 



2ct eq , corresponds to 8, 



<</ 



to 9dec 

0.3°. 

We consider the root-mean-square fluctuation of the temperature of the 
CMB on angular scales 9 ^> 1.7°, i.e. fluctuations that entered the horizon 
after the decoupling of matter and radiation. On these large scales 
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We use a "window function" 
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Figure 8: Histogram of the absolute peculiar velocities of galaxies. 



Then the mean-square fluctuation of mass in randomly chosen windows of 
volume Vw is: 
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(23) 



at early times, i.e. Q ~ 1. For the power spectrum (|20| ) in the range of 
interest of n we obtain approximately 

(24) 
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The fluctuation of the temperature of the CMB radiation is given by 
the Sachs- Wolfe effect P, ||, §|. To obtain an analytic expression we use the 
prescription: 

rmsfi \ 1V1 I rms,Hfi 



T 



where (5M/M) rms H g is the root-mean-square fluctuation of mass on the scale 
corresponding to angle 9 when that mass crossed inside the horizon. The pre- 
scription fl25|) is in agreement with numerical calculations^, ||. From (f25|), 
, (|19D for r s , and replacing A 1 ' 2 by the amplitude at horizon cross- 
A 1 / 2 = a H f{Sl )A 1 / 2 = ^ 2 ^ 1 /(fi )A 1 / 2 , and substituting numerical 



ing a cH 

values, we finally obtain 
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Figure 9: Histogram of the component v y of the absolute peculiar velocity. 



with 6 in radians, and f(Cl ) = 3(1 — Q ~ 1 )/(5A(Q )). The published fluctu- 
ation of the temperature of the CMB is given in terms of the amplitudes of 
spherical harmonis at = (|a™|)[||: 
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for 6 < 2vr. 



3.6 Peculiar velocities 

The absolute peculiar velocities of galaxies with respect to the CMB, shown 
in Figures [8] and |9], have relatively large contributions from wavevectors k p» 
keq/v which are well beyond the range of our simulations. Furthermore, these 
absolute peculiar velocities are difficult to measure. We therefore define a 
local peculiar velocity in spheres of radius r = 9.2h^ x Mpc: 



\ v pec) 1 



1 r 

T7 / j \ v pec,i ~~ \Vpec) 
i=X 



:28i 



where N is the number of galaxies in a sphere and the average is over all 
spheres. This local peculiar velocity can be estimated from the lengths of 
the "Fingers of God".@ See also [§]. 

3.7 The cosmological constant 

Until now we have set ^a — 0. Let us consider a low density spatially flat 
universe with Q = 0.3 and ^a = 0.7. By numerical integration we obtain 
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\Vnp.c/ r 



(ST/T) rmse=01U5 = (1.06 ± 0.04 ± 0.20) ■ 
= (220 ± 50 ± 60)km s" 1 
,~ (748 ± 150 ± 350)km s" 1 
5N/N = 0.84 ±0.17 ±0.10 
r = (5.4±1.0±2.8)/^ 1 Mpc 
7 = 1.77 ±0.04 ±0.54 
ln(0, ■ h 3 ■ Mpc 3 ) = -4.6 ± 0.4 ± 0.5 
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Table 3: Benchmark data[||, [7|, |8| used to define a \ 2 to compare the model 
with observations. The first error is observational, and the second one in- 
cludes theoretical errors, statistical errors of the simulations and errors of the 
fits. We add them in quadrature. All errors in this article are one standard 
deviation. 

H t = 0.964, uo = 0.478, A = 0.0, S m = 1.85 and / = 1.17. The simulation 
then proceeds as before. The results are discussed in the next Section. 



3.8 Constraints on Qq? n an d A 

To compare quantitatively the predictions of the model with observations we 
define a x 2 as indicated in Table |3|. This \ 2 has seven terms and the model has 
three parameters (f2 , n and A), so we are left with four degrees of freedom. 
We have excluded from x 2 the Tully-Fisher and Samurai parameters /i and 
£, which are well satisfied, because they are common to all simulations. The 
theoretical errors of 5N/N, 7 and r$ quoted in Table || are obtained from an 
estimate of the error of v pec in the approximation (|TT1), and the propagation 
of this error as determined from numerical simulations. 

Several simulations are compared with observations in Table f|. For each 
pair (Q , n) we have obtained A from COBE data (see Table |^) and Equation 
(|26|) . Note that we obtain good quantitative agreement with observations for 
several pairs (Q , n). 

The best simulation with Q\ = has \ 2 = 3.2, the critical density 
Qq = 1.0, the scale invariant Harrison-Zel'dovich slope n = 1.00, and A = 
1706Mpc 3 . This simulation has 2053 galaxies, Mi umif = 9- 10 12 M sun n; Mm r2 _1 , 
t> c * = 118km/s, 0* = 0.0065/i^Mpc" 3 , Ri um * = 2.8Mpcfi^ m fi _1 , (v pec ) rms = 
486km/s, 5N/N = 0.744, r = (5.9 ± l.O^Mpc, 7 = 1.5 ± 0.3 (these 
are errors of the fit), and a fraction of mass in galaxies e = 0.97. Several 
distributions of this simulation are presented in Figures to ^. 

Some results of the simulations are shown in Tables || to §. 
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Q , Q A \n 


1.50 


1.25 


1.00 


0.75 


0.50 


0.20, 0.0 


11.8 


11.4 


* 


* 


* 


0.35, 0.0 


72.4 


5.7 


12.8 


* 


* 


0.50, 0.0 


196.6 


5.8 


8.5 


** 


* 


0.70, 0.0 


24.7 


16.2 


6.6 


11.5 


* 


1.00, 0.0 


46.3 


10.7 


3.2 


7.3 


78.0 


1.50, 0.0 


162.7 


98.4 


3.7 


5.3 


14.3 


0.30, 0.7 


38.3 


9.2 


12.2 


* 


* 



Table 4: We compare the model with observations by giving the x 2 f° r 4 
degrees of freedom for several CLo, £l\ and n. The amplitude A 1 / 2 was ob- 
tained from COBE data[0 and Equation (|26"D for each pair (Qo,n). Entries 
with a star are excluded because zero galaxies were generated. Entries with 
X 2 > 9-5 are excluded with 95% confidence. 



ffo, ^a \ n 


1.50 


1.25 


1.00 


0.75 


0.50 


0.20, 0.0 


3017 


708 











0.35, 0.0 


1929 


2584 


397 








0.50, 0.0 


1444 


2229 


1866 


22 





0.70, 0.0 


1146 


1765 


2285 


778 





1.00, 0.0 


890 


1398 


2053 


1937 


194 


1.50, 0.0 


846 


1120 


1635 


2057 


1550 


0.30, 0.7 


1451 


2277 


1094 









Table 5: Number of galaxies in the simulations as a function of Qq, Q\ and 
n. 
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Q , ^a \ n 


1.50 


1.25 


1.00 


0.75 


0.50 


0.20, 0.0 


1.06 


0.46 


* 


* 


* 


0.35, 0.0 


1.31 


0.76 


0.49 


* 


* 


0.50, 0.0 


1.55 


0.93 


0.49 


** 


* 


0.70, 0.0 


1.62 


1.15 


0.64 


0.44 


* 


1.00, 0.0 


1.37 


0.96 


0.74 


0.50 


0.39 


1.50, 0.0 


1.35 


1.40 


0.65 


0.53 


0.38 


0.30, 0.7 


1.34 


0.81 


0.45 


* 


* 



Table 6: 5N/N as a function of Qq, Q\ and n. Entries with a star have zero 
galaxies. 

4 Conclusions 

We have developed a simple, fast and predictive model of the hierarchical 
formation of galaxies. We obtain quantitative agreement with observations 
(within the limitations of the model, i.e. outside of the core of galaxy clus- 
ters). The only free parameters of the model are Qq and the power spectrum 
of density fluctuations, i.e. n and A. The COBE observations determine A 
as a function of Qq and n. The model provides insight into the hierarchical 
formation of galaxies, and is useful to study the onset of galaxy formation, 
the merging of galaxies, the redshift-luminosity distributions, and so forth, 
and to further constrain Q , n, A and N e ff. 

The model is robust for two reasons: i) Galaxies are not stepped forward 
in time so errors do not accumulate in time; and ii) Galaxies are placed 
where the exact solution for the density diverges. Weaknesses of the model 
are: i) Galaxy formation and merging in the model is step-wise rather than 
continuous; and ii) The peculiar velocities and displacements are calculated 
in the linear approximation. Therefore, we expect that the simple model 
will not reproduce the center of galaxy clusters in detail which are in the 
process of merging and have gone non-linear, nor will it predict correctly the 
galaxy-galaxy correlation at small separation. 

Comparison of general galaxy observations with simulations determines, 
with 95% confidence, that fi > 0.25 and 1.1 - 0.3fi < n < 1.4 - 0.2fi 
(assuming Q^ = and N e ff = 3.36). The best simulation has a \ 2 P er 
degree of freedom equal to 3.2/4 and corresponds to Qq = 1.0, ^a = and 
n = 1.00. A low density flat universe with Qq = 0.3 and fl^ = 0.7 is still 
allowed. See Table f| for details. 

In summary, we have developed a simple, fast and powerful tool to study 
the large scale structure of the universe. 
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Q , VL A \n 


1.50 


1.25 


1.00 


0.75 


0.50 


0.20, 0.0 


568 


234 


* 


* 


* 


0.35, 0.0 


916 


428 


204 


* 


* 


0.50, 0.0 


1086 


594 


279 


** 


* 


0.70, 0.0 


1254 


761 


382 


192 


* 


1.00, 0.0 


1531 


967 


486 


284 


130 


1.50, 0.0 


1668 


1085 


645 


397 


220 


0.30, 0.7 


756 


337 


156 


* 


* 



Table 7: (v pec ) in km/s as a function of Q , Q A and n. Entries with a star 
have zero galaxies. 



Q , fi A \ n 


1.50 


1.25 


1.00 


0.75 


0.50 


0.20, 0.0 


43 


33 


* 


* 


* 


0.35, 0.0 


69 


60 


48 


* 


* 


0.50, 0.0 


80 


80 


69 


** 


* 


0.70, 0.0 


97 


97 


85 


77 


* 


1.00, 0.0 


118 


118 


118 


109 


97 


1.50, 0.0 


143 


129 


139 


131 


143 


0.30, 0.7 


67 


58 


53 


* 


* 



Table 8: Velocity of circular orbits of L* galaxies, i> c *, in km/s as a function 
of Qq, fl\ and n. Entries with a star have zero galaxies. 
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